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ABSTRACT 

While limited to low spatial resolution, the next generation low-frequency radio 
interferometers that target 21 cm observations during the era of reionization and prior 
will have instantaneous fields-of-view that are many tens of square degrees on the 
sky. Predictions related to various statistical measurements of the 21 cm brightness 
temperature must then be pursued with numerical simulations of reionization with 
correspondingly large volume box sizes, of order 1000 Mpc on one side. We pursue a 
semi-numerical scheme to simulate the 21 cm signal during and prior to Reionization 
by extending a hybrid approach where simulations are performed by first laying down 
the linear dark matter density field, accounting for the non-linear evolution of the 
density field based on second-order linear perturbation theory as specified by the 
Zel'dovich approximation, and then specifying the location and mass of collapsed dark 
matter halos using the excursion-set formalism. The location of ionizing sources and 
the time evolving distribution of ionization field is also specified using an excursion-set 
algorithm. We account for the brightness temperature evolution through the coupling 
between spin and gas temperature due to collisions, radiative coupling in the presence 
of Lyman-alpha photons and heating of the intergalactic medium, such as due to a 
background of X-ray photons. The hybrid simulation method we present is capable 
of producing the required large volume simulations with adequate resolution in a 
reasonable time so a large number of realizations can be obtained with variations in 
assumptions related to astrophysics and background cosmology that govern the 21 cm 
signal. 
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1 INTRODUCTION 

Observations of the 21 cm line of neutral hydrogen are cur- 
rently considered to be one of the most promising probes of 
the epoch of reionization (EoR) and possibly even the pre- 
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ing to a frequency selection for observations, the 21 
cm data provides a tomographic vi ew of the re i oniza- 
tion and pre-reioniza tion process dSantos et al.l 120051 ; 
iFurlanetto et al]|2004aT ) and can also be an exquisite cos- 



mological probe providing complimentary information when 
compared to cosmic microwave backgro u nd anisotropics 
jMao et all l20Qg|; ISantos fc Cooravl 120061 ; iMcQuinn etafl 
2006; iBowman et all l2007l h New large area radio inter- 
ferometers are now being developed, with the promise of 
measuring this signal, namely, LOFAFlQ in the Netherlands, 
in Australia and the low-frequency extension of the 
planned Square Kilometre Array (SKA0). 

Motivated by the observational possibilities offered by 
the current and upcoming low frequency radio interfer- 
ometers a great deal of effort has been underway in or- 
der to fully understand and generate the expected 21 
cm signal that will be seen by these experiments (see 
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iFurlanetto et al.ll2006l for a review). Analytical models can 
be very useful to quickly generate the signal with high 
dynamic range, in particular the power spectrum of 21 
cm br i ghtne s s temperatu r e fluctuations (IFurlanetto et all 
l2004bl . 120061 ; ISethil 120051 ; iBarkanal 120071 ) and for testing 
the ability of a given experiment to constraint cosmolog- 
ical and astrophysical parameters ([Santos fc Cooravl |2006| ; 
iMcQuinn et al.ll2006l ; iMao et al.ll2008r ). The analytical mod- 
els have also been useful to understand the possible contribu- 
tions to the 21 cm signal at high red shifts |Barkana fc Loea 
120051 ; IPritchard fc Furlanettol 120071 ), although it is at low 
redshifts (z < 10) that they seem to provide a b etter descrip- 
tion of the 21 cm 2-point correlation function (Santos et al.l 
2008). However, these models still have problems dealing 
with the spatial distribution of the reionization process, such 
as bubble overlap and ignores complicated astrophysics dur- 
ing reionization. 

Numerical simulations, on the other hand, can poten- 
tially provide an improved description of the 21 cm bright- 
ness temperature signal from first principles. They typically 
involve a combination of a N-body algorithm for dark mat- 
ter, coupled to a prescription for baryons and star formation 
(often complemented by a hydrodynamical simulation) , plus 
a full r adiative trans f er code to propagate t he ionizing ra- 
diation (lGnedinll2000l; Razoumov et alJ|20oS Sokasian et al.l 
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2007^ ; IMcQuinn et al. 120071 ; IShin et all 120071 ; 
2007l ; lBaek et al.ll2009T ). These simulations have 



the advantage of properly dealing with the spatial distribu- 
tion of the fields (such as bubble overlapping or the non- 
sphericity of the bubbles) and can be used to extract more 
information than just the 21 cm anisotropy power spectrum. 
The disadvantage is that they are slow to run, being con- 
strained in dynamical range to sizes typically smaller than 
143 Mpc, and producing a large set of simulations for cosmo- 
logical and astrophysical parameter estimates is unrealistic. 

Instead the 21 cm community makes use of hy- 
brid ap proaches involv i ng sim ulations and analytical tech- 
niques. I Thomas et al.l |2009T ) generate faster simulations 
by post-processing a dark matter simulation with a 1- 
D radiative transfer code to create bubbles around the 
sources of radiation embedd ed in the dark matter halos. 
iMesinger fc Furlanettol (2007) hav e taken a hybrid approa ch 
between the analytical models (|Furlanetto et al ] l2004bl lah 
and the spatial description provided by the full numerical 
simulations. This is based on 3-D Monte Carlo realizations 
of the dark matter density field combined with an "ana- 
lytical prescri ption" in or d er to generate the halo and ion- 
ization fields. Z ahn et al. (2007) also considered a similar 
analytical prescri ption but based dir ectly on the linear den- 
sity field (see also ! Alvarez et al]|2009l ). Although much faster 
than previous numerica l appr oaches, the hybrid approach of 
IMesinger fc Furlanettol |_2007) is still constrained to "small" 
sizes depending on the amount of shared memory available 
in the computer and to low redshifts since it neglects the 
spin temperature evolution for the 21 cm calculation. 

It is useful to note that the proposed low-frequency in- 
terferometers that plan to observe the 21 cm signal from 
reionization have low spatial resolution capabilities but very 
large fields of view - 5x5 deg 2 or more, requiring boxes of 
at least 1000 Mpc in size if we want to properly simulate a 
single field of view and pass such a simulation through the 



observation pipeline in order to understand how instrumen- 
tal noise and systematic effects impact the observations. At 
the same time we need to be able to resolve halos of the order 
of 10 8 Mq corresponding to a cooling temperature of 10 4 K 
that is deemed necessary to host the ionizing sources respon- 
sible for reionization. This implies a resolution of 0.14 Mpc 
requiring the use of large boxes with (7000) 3 cells, making 
even a hybrid method not useful. Moreover, at high red- 
shifts we need to take into account detailed physics that 
determine the 21 cm brightness temperature, such as the 
heating of the inter galactic medium (IGM) or the Lyman-a 
coupling of the spin temperature to the gas temperature, 
which makes it a challenging exercise to implement a com- 
plete radiative transfer /gas dynamics algorithm in a large 
volume simulation box. 

In this paper, we propose a semi-numerical tech- 
nique, partially following the hybrid approach prescribed in 
IMesinger fc Furlanettol (|2007h , capable of quickly generating 
an end-to-end simulation of the 21 cm signal even at high 
redshifts when the effect of the spin temperature is non- 
negligible. Moreover, this method can be used to simulate 
very large volumes (e.g. 1000 Mpc), crucial to simulate the 
field-of-view of next generation of radio telescopes, without 
sacrificing the speed or requiring unpracticable amounts of 
computer memory. The code to generate this type of simula- 
tion will be provided publicly onlin^B and it will be subject 
to continuous improvement through calibration against full 
radiative transfer/hydrodynamic simulations. The large vol- 
ume simulations created with this method will be part of the 
SKA Design Studies Simulated Skies (S 3 ) initiative and we 
hope this approach can prove useful to generate sky models 
for the future 21 cm experiments, which is crucial to test 
for calibration issues and foreground removal as well as to 
study the features of the 21 cm signal that can be probed 
by a variety of experiments and to check the dependence 
of various statistics extracted from the observations on the 
underlying astrophysical and cosmological parameters. 

In the next Section we introduce the algorithm to gener- 
ate the simulation of the 21 cm signal up to high redshifts, 
while in Section 3 we explain how to extend the simula- 
tions to very large volumes. For results presented here, we 
assumed a flat ACDM universe with the following cosmolog- 
ical parameters: Q m = 0.28, n b = 0.046, h = 0.70, n s = 0.96 
and a s = 0.8 17 based on the resu lts from WMAP5, BAO 
and SN (see iHinshaw et al.l 120091 and references therein) . 
Throughout the paper we quote all quantities in comoving 
units unless otherwise noted. 



2 SEMI-NUMERICAL SIMULATION UP TO 

z=25 

In this Section, we outline our semi-numerical approach, first 
focusing on establishing the linear density field and from it 
the dark matter halo distribution followed by the ionization 
field, with first order perturbations to account for non-linear 
effects. 
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2.1 From linear density fluctuations to the halo 
distribution 

The underlying basis for our simulation is a Monte-Carlo 
realization of the dark matter linear density field assuming 
a Gaussian probability distribution function for the linear 
overdensity as the initial state, which is then later evolved 
and used to find the collapsed structures and the correspond- 
ing ionized bubbles. This i mplementation partially follow s 
the algorithm prescribed in iMesinger fc Furlanettq l|2007t ). 
which is based o n the "peak-patch" approach introduced by 
iBond fc Mversl (| 19961 ). The principle behind this procedure 
is that at the high redshifts such as those where reioniza- 
tion occurs, the spatial distribution of dark matter halos 
can be obtained with linear theory, and non-linear gravita- 
tional effects can be accounted for by using only first-order 
perturbations o n the linear field - i.e., with the Zel'do vich 
approximation l|ZerDovichlll970l ; lEfstathiou et al.lll985t ). 

As part of the simulations presented here, we first start 
by creating a realization of the linear matter density field, us- 
ing a random Gaussian distribution to generate the Fourier 
space density fluctuations in a box with comoving size L and 
N cells: 



5(k) 



Pss(k)L3 



(a k + iby 



(1) 



where at, &k are Gaussian random variables with unit vari- 
ance and zero mean and Pss(k) is the dark matter power 
spectrum. We define the power spectrum P aa (k) for a generic 
quantity a, as (a(k)a(k')*} = (27r) 3 <5f,(k-k')P aa (fc) for con- 
tinuous fields, where So is the Dirac delta function. The 
power spectrum can be determined from linear evolution 
of the primordial density fluctuations. Although numerical 
codes such as CAMB and CMBFAST could be used for a 
more accurate evaluation of the power spectrum, the simula- 
tions p resented in this paper us e the approximation provided 
by the lEisenstein fc Hul (Il998l ) fitting formulae for the dark 
matter transfer function in standard ACDM models (note 
that the our simulation code also includes a routine to read 
in a transfer function from a Boltzmann equation solver such 
as CAMB) 

The density field in real space then follows from the 
inverse Fourier transform: 



<S(x) = ys ^<5(k)expik- x 



(2) 



where x and k are discretized taking into account the simu- 
lation box. This realization of the dark matter density field 
can be used to identify collapsed halos using the excursion- 
set formalism, in which a given region is considered to un- 
dergo gravitational collapse if its mean overdensity reaches 
a certain critical value of order of S c (z) ~ 1.68/D(z), assum- 
ing spherical symmetric collapse (D(z) is the linear growth 
factor normalized so that D(0) = 1). High resolution N - 
body simulations (|Sheth fc Tormenlll999l : ISneth et alj|200lf ) 
showed that by considering ellipticity in non-linear gravi- 
tational collapse, the critical overdensity would become a 
function of the filter scale as: 



S c (M,z) = VaS c (z) 
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where 8 c (z) is the quantity introduced above and a, b and 
c are fitting parameters. We have adopted this definition in 



order to identify dark matter halos. This algorithm simply 
applies a series of top-hat filters of decreasing size to obtain 
the mean density around any given point of the density field. 
Note that in order to make the process faster, we applied the 
filtering procedure as a multiplication in Fourier space (with 
the corresponding Fourier transform of the window function) 
and then transformed back to real space. A given cell is then 
considered to be the center of a halo of mass M when the 
condition 8{M,x) > S C (M, z) is verified. We also consider a 
full exclusion criterion, which means that halos overlapping 
with a previously marked halo are discarded. This algorithm 
allows then to obtain a catalog of halo positions and masses 
which can be used to identify the star formation regions 
responsible for ionizing the surrounding IGM. An important 
factor in this analysis is the resolution used in our linear grid, 
which sets the minimum size of halos that can be identified. 
If one wants to go as far as to resolve halos of mass 1O 8 M0, 
we need to have a resolution of order of 0.14 Mpc. 

The above spatial resolution (which sets the grid size of 
our simulation box), together with the available RAM mem- 
ory of the machine used to create these first simulations, 
justifies the maximum size of the simulation to be of L=300 
Mpc for N=1800 3 cells (resulting in a minimum halo mass 
of 1.65 x 10 8 M Q ). For the results presented here we used a 
machine with 32 CPUs and 64 GB of shared RAM mem- 
ory, although only 20 CPUs were used at any given time. 
This algorithm requires the use of large Fast Fourier Trans- 
forms (we used the fftw-3.1.1 packagqj) to filter the density 
field in Fourier space and thus requires a significant amount 
of allocatable shared memory in order to create simulation 
boxes with enough resolution. In our machine, we managed 
to obtain the full halo catalog in a reasonable running time 
of less than 1 hour per redshift. 

In figure[T]we show the mass functions obtained at z=10 
for two different configurations (L=300 Mpc;N =1800 3 and 
L=143 Mpc;N=1536 3 ) and compa re them with the results 
from the full N-body simulation of iTrac et all (|2008T l which 
covers a region with L=143 Mpc. The strai ght line represents 
the ex pected mass function as derive d bvlSheth fc Tormenl 
( 1999), using the fitting parameters by I Jenkins et al.l l|200lh . 
One can see that both simulations agree quite well with the 
theoretical curve although the number of low mass halos 
in the simulation with the lower resolution (L=300 Mpc) is 
slightly underestimated. However, we have checked that this 
effect has a minor impact on the ionized hydrogen distribu- 
tion and can be seen as a good compromise in order to cover 
the largest possible volume with this algorithm. 

With the halo distribution derived from the linear den- 
sity field, the positions of both the halos and the dark mat- 
ter were corrected to include non-linear dynamics. This was 
done by the means of the Zel'dovich approximation, i.e., 
by using the linear velocity field, obtained from our simu- 
lations, to adjust the positions of both the halos and dark 
mat ter distribution, starting at a n arbitrary high redshift 
( see IMesinger fc Furlanettd 120071 for details on this imple- 
mentation). 

In figure [2] we show the power spectrum at three dif- 
ferent redshifts of the dark matter fluctuations with and 
without the non-linear corrections, in the dimensionless form 
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Figure 1. Mass functions at z=10 taken from our halo filter- 
ing prescription with L=143 Mpc and N=1536 3 (blue circles) 
and L=300 Mpc and N=18 00 3 (black diamo nds). The results 
for the N-body simulation of iTrac et all l|2008l l with L=143 Mpc 
are also s hown as green sq u ares. The line show the mass func- 
tio n from ISheth fc Tormenl \l99$) . using the fitting parameters 
by Ijenkins et al. I l l200lf) . Error bars are included only for the 
L=300Mpc data for viewing purposes. 
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Figure 2. Normalized power spectrum of matter density fluc- 
tuations at rcdshifts 7.5 (upper curves), 10 (middle curves) and 
14 (lower curves). Straight and dashed lines correspond to the 
power spectrum in our simulation, before and after applying the 
Zel'dovich corrections. 

A(fc, z) — k 3 P(k) /2-7T 2 ). In terms of the density power spec- 
trum, we can see that the Zel'dovich approximation in- 
creases the correlations at small scales similar to what is 
seen in full N-body simulations. 

2.2 From the halo distribution to ionization field 

With the corrected halo and density fields, the ionization 
regions can be determined using a similar excursion-set al- 
gorithm. The principle behind this procedure is that the 
galaxies formed inside dark matter halos will produce a given 



amount of photons (dependent of the halo mass) that will 
ionize the surrounding hydrogen generating ionized bubbles. 
The efficiency of this process can be quantified by one pa- 
rameter £, so that a region will be ionized when: 

fcoii > C 1 , (4) 

where f co u is the f raction of mass that has collapsed in ha- 
los in that region l|Furlanetto et al1l2004bl ). We assumed a 
constant efficiency parameter although this can be easily 
extended to include redshift and mass dependence without 
degradation of the speed of the simulation. We note that 
this is an approximate model and contains several limita- 
tions when compared with full radiative transfer methods. 
In particular it does not include radiative feedback effects 
and spatially dependent recombinations. Neverheless, com- 
parison with such more detailed (albeit slower) simulations 
show that our bubble filtering algorithm is a very good ap- 
proximation providing maps with enough accuracy to be 
compared against future 21cm experiments. 

We use this principle to construct our ionization field 
(xi), averaging both the halo mass and density fields with 
the spherical top-hat filters of decreasing size. Whenever the 
condition in equation Q was verified, all cells contained in- 
side the sphere where flagged as completely ionized, allow- 
ing at the same time for overlapping regions. To increase 
the speed of the simulation, we chose to smooth the density 
and halo fields into a 600 3 grid before starting the filtering 
algorithm. The impact of reducing the resolution is to in- 
crease the minimum bubble size, causing some small scale 
ionization to be lost. To account for this effect we considered 
a fractional description of the ionization state within each 
cell. After filtering our box up to the cell size, we searched 
for non ionized cells containing halos, to which we attributed 
an ionization value of Xi = Vi/V ce ii, the ratio between the 
ionized volume which follows directly from our barrier con- 
dition Vi = M co nC/(l + 5i)pm and the cell volume. 

Another important aspect of the degrading resolution 
for ionization calculation can be the overall loss of informa- 
tion on small scales when using smoothed boxes to apply the 
Zel'dovich approximation. The use of small size boxes is de- 
sirable to reduce computational costs, but reducing the box 
size of the velocity field used to displace halos and dark mat- 
ter particles can have a significant impact in the final mean 
ionization fraction (Si). We checked that by performing all 
calculations with smoothed boxes of 300 3 or 600 3 cells the 
value of Xi can change by up to 5% for a fixed value of f . To 
keep the maximum of information at all scales we decided 
to perform all corrections with 1800 3 boxes and smooth 
the density and halo distributions only after applying the 
Zel'dovich corrections. There is no significant increase in 
the computational cost when implementing the Zel'dovich 
approximation on the original box size. 

As a result of this multi-step algorithm, we obtained 
the distribution of ionization bubbles in our simulation box 
which can be used to derive several relevant quantities. In 
figure [3] are shown 4 slices of our ionization field at decreas- 
ing redshifts. These figures exhibit the characteristic bubble 
structure of ionized hydrogen increasing in size at lower red- 
shifts until converging into a completely ionized IGM. The 
four panels of figure [3] correspond to x% =0.08, 0.24, 0.46, 
0.82, with the bulk of the reionization happening between 
z=14 and z~7.5. Obviously, these values are intrinsically de- 
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Figure 3. Slices of the ionization field at four different rcdshifts: clockwise from top-left z = 12, 10, 9, and 

pendent on the assumed efficiency (here we assumed a value 
of 15.1), which could be constrained through observational 
data. This is shown in figure [4] where we compare the evo- 
lution of the mean ionization fraction with redshift using 
( — 6.7 and 15.1. We see that by changing the efficiency pa- 
rameter we can change the redshift where reionization starts 
and fit to a desired optical depth (we can quickly generate 
the ionization field once we change £). 

In order to characterize the statistical distribution of 
the ionization field we computed its power spectrum (P XiXi )- 
In figure [5] we show the results using two configurations of 
L=300 Mpc, N=600 3 and L=143 Mpc, N= 768 3 and com- 
pare it to the results from lTrac et al. (2008) which also has 
L=143 Mpc, N=768 3 . The efficiency parameter was adjusted 
to obtain the same ionization fraction Xi(( — 7 for both 
the L=143 Mpc and L=300 Mpc simulations) . Since these 
configurations correspond to different spatial resolutions, a 
smoothing was performed in the N=768 3 ionization boxes to 
roughly match the resolution of the 0.5 Mpc per cell achieved 
in the 300 Mpc simulation. Two important conclusions can 
be drawn from this figure: the excursion-set formalism al- 
lows to obtain a bubble power spectrum that matches the 
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Figure 4. Evolution of the ionization field Xi using f = 6.7 
(dashed line) and f = 15.1 (solid line). 
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Figure 5. Power spectrum of the ionization field for our sim- 
ulation with L=f43 Mpc,N=768 3 (red dotted curv es), L=300 
Mpc,N =600 3 (solid blue curves) and the simulation of lTrac et al,l 
(2008) (dashed green curves) for 2 = 9 and X{ = 0.15 



one obtained from a full N-body radiative transfer simula- 
tion when the volume size and cell resolution are the same. 
However, an important difference is observed between the 
power spectrum at large scales obtained in our simulation 
with a larger volume (L=300 Mpc) when compared to the 
lower volume ones. We believe this is due to the fact that we 
are finding larger bubbles in the 300 Mpc box, while smaller 
bubbles are absorbed in larger ones (a similar effect was 
observed in lMesinger Sz Furlanettoll2007l ). We have checked 
that such differences do not originate from the lower number 
of low mass halos obtained in the L=300 Mpc simulation or 
the fact that the mass function extends to larger masses in 
the larger simulation. 



2.3 The 21 cm signal 

With ionization maps at different redshifts we can focus on 
the predicted 21 cm signal from neutral hydrogen during 
the pre-EOR and the EOR. In the following Section we will 
briefly present the theoretical aspects behind this radiation 
and the results from our simulation. For a given observa- 
tional frequency, the 21 cm signal corresponds to an inten- 
sity variation along the line of sight of 
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where Ts is the spin temperature of the IGM, z is the red- 
shift corresponding to the frequency of observation (l + z = 
v it /v), where v 2 i =1420 MHz) and T 7 =2.73(l + z)K is the 
CMB temperature at redshift z. The optical dep_th corre- 
sponding to the hyperfine transition is given by (|Fieldll 19591) : 

3c 3 hA 10 n H i 

T ~ Wku2iTs(l + z)(dV r /dr) [ ' 

where A\o is the spontaneous emission coefficient for the 
transition (2.85 x 10~ 15 s" 1 ), nm is the neutral hydrogen 
number density and dV r /dr is the comoving gradient of the 
total radial velocity along the line of sight. Neglecting pe- 
culiar velocities, one has dV r /dr = H(z)/(l + z). The neu- 



tral density can be expressed as nm = fniXpb/m p where 
Jhi = Phi/ph is the fraction of neutral hydrogen (mass 
weighted), X « 0.76 is the hydrogen mass fraction, pt is 
the baryon density and m p the proton mass. Since the bary- 
onic content follows the dark matter distribution, the 21 cm 
temperature can then be written as: 



8T b (v)& 23 fm(l + S) 1 



h_ 

07 



1 'nth 2 



0.15 



l + z 
10 



1/2 



1 + l/Hdv r /dr 



0.02 
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where dv r /dr is the comoving gradient of the line of sight 
component of the comoving velocity (obtained directly from 
the density field using a fourier transform). In high density 
regions of the simulation one can have 1/H dv r /dr — 1 and 
therefore we imposed a limit to this gradient of dv r /dr > 
—0.7H in order to avoid singularities. Note that this only 
affects a very small number of cells (< 0.05%) and we verified 
that in neutral regions we have \dv/dr\ < H . In the above 
equation, both fm and S (as well as the velocity gradient) 
can be obtained as a direct result of our simulation (we used 
/hi = 1 — Xi), but the term depending on Ts is sensitive 
to variations of the spin temperature value over the volume. 
For z < 10, one can make the reasonable assumption that 
T s ^> T 7 since the IGM gas temperature should be well 
above 100K at these redshifts but at higher redshifts the 
contribution to the signal from the spin temperature has to 
be taken into account which we address further below. 

In figure [6] we plot the brightness temperature power 
spectrum for two different stages of the EoR (xi = 0.24 and 
Xi — 0.81). In each figure we also plot the different contri- 
butions to the total power spectrum. This was done by us- 
ing exclusively density, ionization and velocity fluctuations 
at each time in eq[7] One can observe that at the initial 
stages of reionization (low Xi), the main contributions to 
the 21 cm signal come from the density and ionization fluc- 
tuations, meaning that the hydrogen spatial distribution is 
moderately important due to the small ionization bubble 
size. At later stages, ionization bubbles have grown to a size 
where they completely dominate all the other contributions. 
As for the velocity fluctuations, they are sub-dominant at all 
stages, although being more important for higher redshifts. 



2.4 Spin temperature at high redshifts 

At high redshifts (z > 10) the spin temperature (Ts) is no 
longer high enough to saturate the effect in equation [7] and 
we need to take into account the contribution of fluctua- 
tions from the spin temperature to the 21 cm signal. These 
originate from fluctuations in the coupling between the spin 
temperature and the gas temperature and the perturbations 
in the gas temperature itself (Tk) and we can write: 



^ T^ %tot 

Ts 1 + xtot 



Zz. 



(8) 



where x tot = x a + x c is the sum of the radiative 
and collisional coupling parameters. Collisions can be 
important for decoupling the HI 21 cm spin tempera- 
ture from the CMB, e specially at very high redshifts 
(z > 30) E 3 d2005h and are straightforward to ap- 
ply to the simulation jAllison fc Dalgarnol 1 19691 : IZygelmar] 
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Figure 6. Brightness temperature power spectrum (solid blue line) and its different contributions at redshifts z=10 and z=8 with 
£ = 15.1. The green dashed-dotted lines correspond to the power spectrum taking only velocity fluctuations, the yellow dashed lines only 
have density fluctuations and the red dotted lines only for ionization fluctuations. 



20051; iFurlanetto fc Furlanettd 120071 ; iKuhlen et all [2006; 
Hirata fc Sigurdsonl 120071 ). The radiative coupling due to 



the absorption of Ly a photon s (the Wouthysen-Field effect, 
I WouthuvserJ 19521 ; lFieldlll959h ■ on the other hand, should be 
dominant for z < 25 and we shall concentrate on calculating 
this effect here. Note however that if we use a model where 
the Ly a coupling is smaller, then the collisional coupling at 
lower redshifts can beco me important close to the sources 
due t o the X-ray emission IjZaroubi et al . 20071: iKuhlen et al.l 
2006). The radiative coupling is given by 

' / " (9) 



Jc 



with J c « 5.552 x 10" 8 (1 + z 



_1 Hz J sr 1 and S a is 



a cor rect ion fact or of ord er u nity jChen fc Miralda-Escudel 



2004 iHiratal l20od IChuzhov fc Shapiro! 



Furlanetto fc Pritchardl 12009 ; IFurlanetto et alj [200' 



We follow the prescription in ISantos et al. 



to 



calculate J a (the spherical average of the number of Lya 
photons hitting a gas element per unit proper area per unit 
time per unit frequency per steradian), given by a sum over 
the hydrogen levels n, 



J Q (X, z) 



(! + *)' 
47T 

4n 



^2 f™c( n ) 



(10) 



*(n) 



dx ip(x + x.',z) e a (v ) 



where we used n max = 16, frcc(n) is the fraction of Lyman- 
n photon that cascade through Lya and the redshift z' in 
equation 1101 is such that x = J z cH~ 1 dz". The emitted 
frequency at z 1 is given by 

- 2 , (! + *') 



(1 + *) 



(11) 



in terms of the Lyman limit frequency i/ll and x max (n) cor- 
responds to the comoving distance between z and z ma x(fi) 
such that: 



l + z n 



(1+z 



[l-(n+ I)" 2 ] 



(12) 



(1-n- 2 ) 

In equation 1101 %p(x,z) is the comoving star formation 




Figure 7. Average comoving star formation rate as a function of 
redshift. Black solid line is from our simula tion and re d dash ed- 
dotted line is from an N-body simulation in IShin et al.1 d2007t) . 



rate density and t a {y) gives the spectral distribution func- 
tion of the sources (defined as the number of Ly n photons 
per unit frequency emitted at v per baryon in stars). We 
can easily assume any model for e a (y) in our code to test for 
different sources of radiation. For this simulation we used: 
Av~ a , with a = 0.9 and A set such that we get a total 
of 20000 Lya photons per bary on which give s a Ly n emis- 
sion similar to the one used in ISantos et al.1 (|200ct ). Note 
however that the amplitude is degenerate with the normal- 
ization of the star formation rate. In order to obtain the 
star formation rate density from the simulation we used the 
halo catalogues generated to create boxes of collapsed mass 
(M co u(x., z)) for two closely separated redshifts and calcu- 
lated: V( x i z ) = f* , §iM B oii(x,z), where /* was used to nor- 
malize the star f ormation rate density to the one used in 
IShin et al.1 (|2007T l (see fig. [TJ . 

Finally, to make the code faster and taking into account 
that we only have the star formation rate density boxes for 
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a range of redshifts, we divided the integral over dQ! dx' 
above in eq. [10] into redshift bins, e.g. 



as well as how much is used to heat the IGM, by integrating: 



i — imax(n) 



E 



d 3 x'Ki(x') ipi(x. + x.'), 



(13) 



with kernel Ki(x') = e a (i/)/(47T x' 2 ) if Zi ^ z < z<+i corre- 
sponding to the redshift bin for which we assume the ipi box 
to be constant and Ki(x') = otherwise. Here, imax(n) is 
set by x ma x(n) but to avoid having to sum over too many 
boxes we made the assumption that when x > 71 Mpc the 
star formation rate was homogeneous. We have checked by 
changing this limit that it is a good approximation (note 
that the star formation rate density fluctuations are quite 
small on large scales, k < 0.1 Mpc -1 ). Moreover, since this is 
a convolution, the expression above can be easily calculated 
with resort to a Fast Fourier Transform (FFT): 



FFT (FFT[Ki\ 

i=0 



FFT[tpi] 



(14) 



FigureJS] left panel, shows the evolution of the radiative cou- 
pling with redshift which is already larger than 1 for z < 18. 
The right panel on the other hand shows the power spec- 
trum of the quantity x a /(l + x a ) which is what contributes 
to the 21 cm signal (we divided the power spectra by the 
square of the average for later comparison, although this 
does not change the amplitude much). We see here that the 
fluctuations peak around z — 20 when x a « 0.1. 

The gas in the IGM is also heated as the reionization 
progresses. In our simulation, we assume that this heating is 
done by x-rays with an emission c onnected to the star for- 
mation rate (see lSantos et al]|2008l ). Our starting point was 
to calculate the total x-ray energy per unit time deposited 
in a given cell: 



exn(x, z) = (1 + z) 2 y — X 

i 

d 3 x'tp(x + x, z') K(x'), 



where 



K(z') 



dv -pr^L ex He t{z,x ) (hv - hv\b) 
47J- x' |z 



(15) 



(16) 



and rii (i =HI, Hel, Hell) is the number density (due to our 
prescription for ionization, we assumed that fnei = Jhi 
and there was no Helll at these high redshifts). The opti- 
cal depth r should in fact depend on the path between x' 
and x but to increase the speed of the calculation we as- 
sumed m to be homogeneous in the integral used to obtain 
the optical depth. Moreover, we tabulated the kernel values 
(which only depend on x') and interpolated them when do- 
ing the FFT making the code faster by a large factor. In the 
equation above, tx(y) is again assumed to be a power law, 
with a spectral index and amplitude compatible with what 
is observed for starburst galaxies (note that this amplitude 
is also degenerate with the star formation rate normaliza- 
tion). We then calculate how much of this energy is used for 
secondary ionizations in the neutral IGM, 



da; e 



= ex. 



/ion 
E t h 



dT K _ 2Tk dn 2e x 
dt 3n At 3fes 



/he 



(18) 



(17) 



where /i on is the fraction of energy converted into 
ionizations and /heat the fraction converted into heat 
( Shull fc van Steenb erg 1985) which depends on the fraction 
of free electrons, x e . We started the integration at z = 25 
assuming an initial temperature of T = 14K consistent with 
adiabatic cooling of the gas from z = 150 (where it is ap- 
proximately equal to the CMB temperature) . Figure [9] left 
panel, shows the evolution of the gas temperature with red- 
shift, where we can see that most of the IGM is heated 
above 100K for z < 11. The spin temperature starts very 
close to the CMB one (no signal) and approaches the gas 
temperature as x a increases, following it for z < 17. In the 
right panel we show the power spectrum of the quantity 
1 — Tcmb/Tk which contributes directly to the fluctuations 
in the brightness temperature (but in this case we divided 
by the square of the average) . Comparing to the right panel 
in fig. [8] we see that fluctuation in x a are dominant at very 
high redshifts (z > 18 when x a < 1) while the perturba- 
tions in the gas temperature are more important at lower 
redshifts, peaking at z — 14 if we multiply the fluctuations 
by the average of 1 — Tcmb /Tk which is what shows up in 
the 21 cm signal. 

Finally, putting it all together we can calculate the 21 
cm signal using equation [7] Calculation of x a and the gas 
temperature takes around 5 minutes for each redshift in our 
machine using 20 CPUs. This obviously depends on how 
many star formation rate density boxes we use for the inte- 
gration which in turn depends on the assumed redshift bin 
(in our case we used dz = 0.25 for each box). In figure [9] left 
panel, we can also see the evolution of the average brightness 
temperature with redshift, which is observed in absorption 
down to z = 12 where it becomes approximately zero when 
Tk ~ Tcmb- 

In figure [TU] we show maps of the 21 cm signal for a 
few very high redshifts taking into account all the effects. 
At high redshifts (z = 21) most of the temperature is zero 
(T s ~ Tcmb) but we already start to see some cold spots 
where the Lya sources couple the spin temperature to the 
cold gas temperature. At z = 17 the spin temperature is 
basically following the gas temperature everywhere and we 
can see fluctuations in the 21 cm signal due essentially to 
the gas temperature (although most of the Universe is still 
cold at this epoch). At z = 14 we can already see the higher 
temperature regions surrounding the light sources located 
in the halos. Finally at z = 11 basically all the IGM is 
heated above the CMB and we start seeing the ionization 
bubbles (xi w 0.1, but note that some of the dark regions 
here are not due to the ionized bubbles but because Ts = 
Tk ~ Tcmb)- Figure mi shows the corresponding power 
spectrum. Note again that at z » 11 the fluctuations in 
the gas temperature still make an important contribution 
to the 21 cm signal, especially on large scales (increasing 
the power), even though the IGM temperature is already 
above the CMB one (Tk ~ 150K and T C mb = 32. 7K). 
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Figure 8. Left: Evolution of the < x a > with redshift. Right: power spectra of x a /(l + x a ) divided by the square of the average of this 
quantity. 
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Figure 9. Left: The evolution of the gas temperature, spin temperature and brightness temperature with redshift. Solid black line - 
CMB; Red dashed - gas temperature; Blue dotted - spin temperature; Dark blue solid/dashed (below) - brightness temperature. Right: 
power spectra of the factor 1 — Tqmb/Tk divided by the square of its average. 



3 EXTENDING THE SIMULATION TO VERY 
LARGE VOLUMES 



Although quite fast and based on first principle physics, 
the semi-numerical simulation presented so far still imposes 
some limitations regarding the volume it can cover. In the 
absence of large computing resources, it is not feasible to run 
a simulation covering a large dynamical range. For exam- 
ple the simulation presented above, covering 300 3 Mpc 3 , re- 
quired to push the limits of our machine, using 20-processors 
and almost 64 GB of shared RAM in order to perform the 
FFT filtering. The main problem with the algorithm that 
was presented is that it still requires a very large number of 
cells in order to achieve sufficient resolution to filter scales 
corresponding to 1O 8 M0 halos and cover large volumes at 
the same time. If we want to simulate the very large volumes 



that will be surveyed by the next generation radio telescopes 
a less demanding algorithm is required. 

This motivated us to implement a slightly simpler 
method when identifying the halos. Using the same resolu- 
tion of N=1800 3 , but with L=1000 Mpc (which corresponds 
to roughly a 5 x 5 deg 2 sky at z=15) the minimum cell size 
of L/N ~0.51 Mpc corresponds to a minimum halo mass of 
~ 1 x 10 10 Mq . Clearly, this configuration does not allow us 
to resolve halos down to the required minimum mass of 10 8 
Mq. However, one can use the expected halo mass function 
to place these smaller halos in our large volume sim ulation 
cells by using the prescription of lWilman et al.l (|2008l ). It fol- 
lows directly from the definition of the mass function that 
the quantity n = dn/dM AM AV corresponds to the mean 
number of halos found in a given cell with comoving vol- 
ume AV. If one would simply perform a Poisson sampling 
for each cell labelled (i,j,k), with mean n(i,j,k), the distribu- 
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Figure 10. Maps of the 21 cm brightness temperature at very high redshifts from the simulation, including all sources responsible for 
brightness temperature fluctuations. 
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Figure 11. The power spectrum of the 21-cm signal at high 
redshifts including all fluctuations in the spin temperature. 



tion of obtained halos would yield a correct mass function 
but their positions would be completely random. In order 
to correlate the halo positions with the underlying density 
field, a normalized bias term can be added, 

n = Ke b{z ' M) ^ z) ^-AMAV , (19) 

where b(z, M) is the bias model for halos at a given red- 
shift and S(z, M) the density field at the same redshift. The 
constant K is a normalization constant that ensures the con- 
sistency of the above expression with the mass function 
when averaging over the large volume box. For a Gaussian 
field, this requirement leads to: 

K _ g-bO.AOV*) 2 /^ (20) 

Here, a(z) 2 is the variance of the density field at a typical 
scale of cell size. The above equations introduce a depen- 
dence of the halo locations on the density field, by amplify- 
ing the overdense regions with respect to underdense ones. 
As for the bias term, it describes the difference in cluster- 
ing between dark matter halos and the mass density field. 
We consider a linear bias developed from the Press- Schechter 
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M (Mo) 

Figure 12. Mass function at z=10 for our simulation with 
L=300Mpc and N=1800 3 ( green squares) comparing with the 
large volume simulation with L=1000 Mpc and N=1800 3 (blue 
diamonds). Halos with masses M < 1O 1O M0 were placed using 
the Poisson sampling method introduced in this section. 

forma lism bv lMo fc White! (| 19961 ) and later modified bv ljind 
(|l99St ): 

/ _ 1 \ / 1 \ 0.06-0.02n 

^M)=(l + ^-)(^- I+ l) , (21) 

where v = 8 c /u(M) and n is the spectral index of the pri- 
mordial density field fluctuation power spectrum. 

We have then used this prescription in complement with 
the excursion-set formalism to perform a dark matter sim- 
ulation with L=1000 Mpc. Since the cell size defines the 
lower mass scale using the exclusion set formalism, halos 
with mass above 10 10 M were obtained using this method. 
For scales bellow the cell size, we did a Poisson sampling 
of the entire box using eq. [T5] for the distribution mean. In 
this way, we allow for the same cell to contain more than 
one halo, as long as the sum of their associated volumes 
does not exceed the cell volume AV. In figure 1121 we plot 
(blue diamonds) the obtained mass function for this run 
and compare it to the higher resolution, L=300 Mpc simu- 
lation obtained using only the excursion set formalism. We 
can clearly see that only this method allows to have a high 
compatibility with the expected mass function with low res- 
olution, whereas the full excursion-set formalism only has 
similar results when the cell resolution is considerably high. 

Despite the good agreement with the mass function, a 
more important aspect is to see how well this formalism can 
reproduce the spatial distribution of the halos and correla- 
tion with the density field. This is a crucial factor, since a 
correct determination of the reionization bubbles depends 
on the relation between these quantities. One way of testing 
for this is by computing the halo power spectrum, P(k)^, of 
the quantity Shh = M co u/M co u — 1. In figure [T3l we show the 
obtained power s pectrum, an d comp are it with the one from 
the simulation of iTrac et al.1 (|20od ). The halo power spec- 
trum derived from our method shows the same characteristic 
shape as the one obtained with the full N-body simulation. 
Once we introduce the first order Zel'dovich corrections for 
the halo positions (red dashed curve) the power spectrum 
becomes almost in perfect agreement with the simulation 
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Figure 13. Halo power spectrum of the halo field obtained 
in our very large volume simulation before and after applying 
the Zel'dovich corrections (dotted dashed and dashed curves, re- 
spective^, compared with the result from the L=143 Mpc of 
ITrac et al'l fcOQgf) . 

(we also checked that this agreement is true for different 
mass bins). Such result shows that this somehow ad hoc 
method can retrieve the "correct" statistical properties of 
dark matter halos. 

Encouraged by this result, we proceed to derive the ion- 
ization bubbles using the same excursion-set formalism as in 
Section 2. For this, we have also applied the Zel'dovich ap- 
proximation for the density and halo field boxes. In order 
to keep a moderate computational time, we used again a 
N=600 3 box for determining the ionization field, which re- 
sulted in a cell size of about 1.5 Mpc. This required the 
need to identify smaller ionized regions (~0.5 Mpc) which 
was done using the same technique introduced in Section 2, 
with the halo catalog being used to identify the ionized vol- 
ume within each non-ionized cell. The results for the power 
spectrum of the ionization field are shown in figure [14] com- 
pared to those obtained with the full excursion-set method 
with L=300 Mpc and for two different efficiency parameters 
at z=9. The curves seem to agree reasonable well showing 
convergence with the 300 Mpc box, although if we go to 
even larger ionization fractions we still find slightly larger 
bubbles in the lGpc box (the largest bubble for the lGpc 
box with Xi = 0.55 has a size of 20 Mpc). Note that the 
bubble size will also b e affected by absorption systems (e.g. 
I Alvarez fc Abell2010l ). 

To illustrate the application of this method we show 
the final brightness temperature maps for the expected 21- 
cm signal in the L=1000 Mpc simulation (fig. [15}. Note that 
we are considering partially ionized cells in this simulation, 
so even the large patches of emission in the figure have a rea- 
sonable amount of ionizing sources, although this is difficult 
to see given the small size of each "pixel" . 



4 SUMMARY 

We presented a semi-numerical method capable of quickly 
generating end-to-end simulations of the 21 cm signal even 
at the high redshifts where the spin temperature is non- 
negligible. The algorithm allows to generate brightness tem- 
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MMpch- 1 ) k(Mpch->) 

Figure 14. Power spectrum of the ionization field for the very large volume simulation (green dashed curves), compared to the same 
result for the L=300 Mpc simulation. Both plots are for z=9 but using different efficiency parameters =5.0 and £=14.0) so that xi = 0.1 
(left) and Xi = 0.55 (right), respectively. 
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Figure 15. The temperature maps for the L=1000 Mpc simulation for two different ionization states at z=9 (as in figure ITU . 



perature boxes with very large volumes, e.g. (1000 Mpc) 3 , 
crucial to properly simulate the field of view of the next gen- 
eration of radio-telescopes, without sacrificing the speed or 
requiring unfeasible computer resource^. 

The corresponding code (SimFast21) implements the 
following prescription: 

• Monte Carlo generation of the dark matter density and 
velocity fields at z=0 assuming Gaussian distribution func- 
tions. 

• Determination of the halo catalog using an excursion- 
set formalism on the linear density field at any given red- 
shift. For large volumes a biased Poisson sampling is used 
to introduce small mass halos. 

6 During the review process of this paper another algorithm was 
presented bv lMesinger et all j201fjl . 



• Adjustment of the halo and dark matter locations using 
the Zel'dovich approximation. 

• Simulation of Ki(x, z) from FFTs of the previous boxes. 

• Extraction of the star formation rate density from the 
halo mass distribution. 

• Extraction of gas temperature fluctuations Tk(h,z) . 

• Calculation of the Ly a coupling and collisional coupling 
parameters. 

• Determination of the 21 cm brightness temperature 
boxes using the above quantities. 

The computational time required for the above steps de- 
pends mainly on the volume to be covered (our bubble filter- 
ing procedure is now piratically independent of the ionization 
fraction due to the FFT based algorithm we are now using). 
This basically affects the halo determination running time, 
since larger volumes imply a larger number of halos and op- 
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orations for overlap checking. Nevertheless, the typical time 
for the computation of the above steps is about 2 hours for 
the simulation with L=1000 Mpc, which can be considered 
remarkably fast. We are expecting that further optimization 
of the code can still reduce this time. 

Although much faster than hydrodynamical numerical 
simulations, our analysis shows that relevant quantities such 
as the halo mass function, the halo mass power spectrum, 
the star formation rate or the ionization fraction power spec- 
trum are all consistent with the numerical simulations with 
which we compared our results to. In particular, the ioniza- 
tion fraction power spectrum is similar to the one obtained 
with radiative transfer codes (RT) when comparing to the 
same volume. We still find some differences on the largest 
scales when comparing to simulations with relatively small 
volumes (143 Mpc), which we can relate to the production of 
larger HII bubbles. These differences become much smaller 
as we increase the simulation volume and the change of size 
in the largest bubbles becomes less important when com- 
pared to the overall size of the simulation. 

The dependence on the astrophysical parameters of the 
simulation was encoded in three functions: the ionization 
efficiency, (, the Ly a spectral distribution function of the 
sources, t a and the X-ray spectral distribution function, ex ■ 
These functions can be easily changed for a model of our 
choice and the code can then quickly generate new simula- 
tions of the signal (even faster if we keep the same cosmol- 
ogy) . By combining all the unknowns into a physically mean- 
ingful small set of parameters we can easily probe the huge 
intrinsic parameter space available at high redshifts probed 
by 21 cm observations. This versatility is included in our 
algorithm. The possibility to generate reionization simula- 
tions from scratch without the need for supercomputers is a 
great advantage allowing to experiment with different mod- 
els. Moreover, this allows room for improvement through 
calibration with full numerical simulations. The code can be 
a useful tool to generate sky models for future 21 cm experi- 
ments, important to test the observation pipeline, study how 
the foregrounds affect the observations and can be separated 
given a noise model, and to develop optimal estimators for 
signal extraction. The code will be released for public use 
in due course and we welcome the community participation 
in its updating and upgrading as well as development of ad- 
ditional applications on observational probes of reionization 
beyond the 21 cm observations. 
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